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This  MEANS  initiative  research  project  consists  of  a  collaborative  effort  between 
students,  postdocs,  and  faculty  at  the  University  of  Illinois  and  at  the  University  of 
Michigan.  The  work  involves  the  integration  of  length  and  time  scale  spanning 
computational  methods  of  investigation  into  a  suite  of  design  tool  for  materials 
optimization.  Concurrent  experimental  measurements  serve  to  motivate  and  validate  the 
simulation  approaches.  Specifically,  the  challenges  posed  by  the  design  of  autonomously 
healing  polymer  matrix  composites1  were  chosen  as  the  test  case  for  the  development  of 
this  computational  framework,  which  includes  reactive  molecular  dynamics  simulations, 
coarse-grained  particle-based  simulations,  and  cohesive-volume  finite  element 
calculations.  The  integration  of  the  different  methodologies  is  achieved  through  data 
exchange  and  output  overlap  matching.  The  project  includes  further  development  of 
specific  numerical  techniques,  implementation  of  an  efficient  coupling  mechanism,  and 
extensive  testing  and  validation  through  experiments. 

1.  Objectives 

The  purpose  of  this  project  is  to  create  a  multi-scale  computational  framework  for  the 
design  of  complex  materials  systems  by  integrating  molecular  dynamics  simulations 
based  on  a  reactive  force  field,  cohesive/volumetric  finite  element  calculations,  and  a 
particle-based  coarse-grained  model  to  bridge  the  gap  between  atomistic  and  continuum 
modeling  techniques.  The  particular  challenge  in  developing  this  computational 
framework  arises  from  the  necessity  to  simulate  the  materials  response  to  mechanical 
loading  on  a  macroscopic  timescale  while  the  structural  developments  associated  with  the 
healing  processes  require  resolution  of  mechanisms  at  atomic  time  and  length  scales. 

While  the  focus  for  this  project  is  on  self-healing  materials  systems,  the  fundamental 
research  pursued  in  this  project  will  yield  tools  of  analysis  and  experimental  databases 
that  impact  materials  applications  beyond  self-healing.  The  multi-scale  simulations  that 
we  are  developing  have  utility  in  the  design  of  composites  and  multiphase  materials  for 
enhanced  fatigue  life  and  increased  fracture  resistance.  The  molecular  dynamics 
simulations  of  in  situ  reaction  and  transport  phenomena,  while  examining  the  stiffness 
and  strength  of  the  resulting  structures,  have  broad  applicability  in  the  field  of  polymer 
processing  and  process  modeling,  especially  in  the  development  of  structural  materials. 

2.  Approach  and  Methodology 

The  design  of  self-healing  materials  systems  requires  very  detailed  and  careful 
optimization  of  processes,  including  the  choice  of  materials  and  repair  agents  based  on 
their  physical,  chemical,  and  kinetic  properties,  as  well  as  geometric  and  structural 
considerations  at  various  levels.  On  the  one  hand,  the  properties  of  a  material  are  very 
fundamentally  rooted  in  the  structure  and  dynamics  of  its  atomic-scale  constituents.  On 
the  other  hand,  the  actual  materials  performance  depends  on  a  complex  scale-spanning 
hierarchy  of  such  elementary  mechanisms.  Incorporating  atomically  resolved  processes 
in  a  modeling  framework  therefore  provides  for  greater  realism  and  accuracy  in 
predicting  materials  behavior.  At  the  same  time  one  needs  to  keep  track  of  how 
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progressions  of  elementary  processes  compound  into  macroscopic  observables. 
Computationally,  this  requires  one  not  only  to  encompass  the  pertinent  length  scales  and 
spatial  resolution,  but  importantly,  the  time  scales  corresponding  to  the  reaction  between 
molecular  building  blocks  and  to  the  fatigue  crack  propagation  must  be  reconciled. 

To  create  a  suitable  multi-scale  design  tool  we  are  combining  explicit-atom  molecular 
dynamics  (MD)  simulations  with  continuum-scale  cohesive/volumetric  finite  element 
(CVFE)  calculations.  Furthermore,  to  overcome  the  time-scale  disparity  between  these 
two  methods,  we  have  developed  a  coarse-grained  particle-based  (CGP)  simulation 
procedure  that  allows  one  to  simulate  atomic-scale  processes  in  an  accelerated  fashion  by 
abstracting  structural  detail  at  the  building  block  level.  CVFE  is  used  to  model  fatigue 
crack  propagation  at  the  macro-scale.  This  part  of  the  project  is  carried  out  in  the  group 
of  Prof.  Geubelle  at  the  University  of  Illinois,  while  our  group  at  the  University  of 
Michigan  is  working  on  the  atomistic  and  intermediate-scale  simulations.  MD 
simulations  are  used  to  simulate  the  structural  aggregation  upon  curing  of  the  healing 
agent,  determine  reaction  rates,  and  establish  the  relationship  between  mechanical 
properties  and  molecular  structures  at  various  stages  of  the  curing  process.  Molecular- 
scale  simulations  furnish  the  loading  rate  dependent  constitutive  behavior  need  as  input 
for  the  cohesive  volume  elements  continuum  calculations.  Conversely,  CVFE  establishes 
the  local  conditions  under  which  the  molecular  state  of  the  material  is  to  be  evaluated. 

While  it  is  generally  possible  to  simulate  molecular  structures  large  enough  to  yield  local 
properties  that  are  representative  of  the  macroscopic  materials  behavior,  a  serious  limita¬ 
tion  in  connecting  atomistic  and  continuum  scales  arises  when  structural  changes  result 
from  macroscopic  mechanical  loading  that  is  many  orders  of  magnitude  slower  than 
atomic  motions.  Particularly,  the  macroscopic  failure  process,  curing  rate  and  deforma¬ 
tion  rate  in  self-healing  composites  compete  on  time  scales  not  attainable  by  MD  simula¬ 
tions.  Moreover,  the  deformation  and  failure  may  induce  structural  defects  whose  spatial 
extents  exceed  those  of  MD  simulation  boxes.  Therefore,  we  insert  an  intermediate-scale 
lattice-based  modeling  approach  to  bridge  between  atomistic  and  continuum  scales. 

The  development  of  our  simulation  framework  is  closely  coupled  with  the  experimental 
work  of  Profs.  White  and  Sottos  at  the  University  of  Illinois.  This  collaboration  allows 
us  to  identify  the  design  issues  that  need  to  be  addressed  by  our  numerical  tools  and 
provides  us  with  experimental  data  to  properly  validate  our  modeling  efforts. 

3.  Accomplishments 

Our  molecular-scale  simulation  approach  encompasses  three  levels:  (i)  a  fully  atomistic 
level,  where  all  atoms  contained  in  the  system  are  explicitly  represented.  At  this  level  we 
determine  reaction  rates  and  structural  information,  (ii)  a  first  coarse-graining  level,  in 
which  larger  polymeric  networks  are  generated  and  the  relationship  between  structure 
and  mechanical  properties  are  determined,  and  (iii)  a  second  coarse-graining  level,  where 
reaction  rates  and  deformation  rates  are  reconciled  and  the  overall  mechanical  response 
of  the  healing  agent  under  application  conditions  is  determined.  The  development  of 
individual  simulation  techniques,  i.e.,  the  reactive  MD  simulations  and  the  CGP  model, 


has  been  completed.  The  methods  have  been  applied  to  study  the  autonomous  repair 
behavior  of  self-healing  polymer  composites  under  strain. 

3.1  Progress  in  Atomic-Scale  Simulations 

Polymerization  and  cross-linking  reactions  have  been  successfully  reproduced  in  atomic 
scale  simulations  for  systems  of  several  hundred  monomers.  These  simulations  revealed 
the  reaction  mechanism  and  polymerization  rates.  The  dicyclopentadiene  (DCPD) 
molecule  posses  two  reactive  sites,  i.e.,  the  C=C  double  bonds  on  the  buckled  hexagonal 
ring  and  the  pendant  cyclopentenyl  ring  (see  fig.  1).  Both  rings  are  subject  to  ring¬ 
opening  methathesis,  which  is  invoked  by  a  catalyst,  e.g.,  the  transition  metal  complex 
Ru(arene)Cl2(PCy3).2’3  The  MD  simulations  using  our  reactive  force  field,4’5  reproduce 
both  ring  opening  mechanisms.  Our  simulations  show  that  because  the  hexagonal  ring  is 
subject  to  significant  strain,  it  is  considerably  more  reactive  and  hence,  the 
polymerization  initially  leads  to  linear  chains.  Cross-linking  sets  in  at  a  slower  rate,  as 
controlled  by  the  opening  of  the  cyclopentenyl  rings.  This  confirms  the  prevalent  views 
concerning  this  mechanism  as  described  in  the  literature.2  Importantly,  through 
simulation  of  the  detailed  bond  formation  processes  we  are  able  to  generate  a  realistic 
model  of  a  polymeric  DCPD  network,  and  establish  its  structure  property  relationships. 
In  particular,  using  this  level  of  simulation  we  can  determine  the  reaction  rate  coefficient 
that  describes  the  polymerization  reaction.  This  only  requires  a  few  hundred  monomers, 
for  as  long  as  all  steric  arrangements  that  can  occur  during  the  polymerization  are 
represented  in  the  simulations.  This  is  readily  achieved  by  studying  a  series  of  structures 
with  different  initial  configurations. 

Explicitly  accounting  for  hy¬ 
drogen  in  the  simulations  is 
costly  in  terms  of  computing 
times,  not  only  because  of  the 
mere  number  of  atoms  to  keep 
track  of,  but  because  of  the 
small  mass  of  hydrogen,  the 
integration  time  step  needs  to 
be  drastically  reduced  for  the 
simulation  scheme  to  produce 
accurate  trajectories.  Struc¬ 
turally,  hydrogen  plays  a  sec¬ 
ondary  role,  e.g.,  it  has  negli¬ 
gible  load  bearing  function, 
and  it  is  therefore  common 
practice  to  account  for  hydro¬ 
gen  atoms  implicitly  by  recon¬ 
figuring  CH  and  CH2  groups 
into  unified  particles,  as  shown  in  fig.  1.  These  new  units  are  characterized  by  sizes  and 
chemical  interactions  representative  of  the  steric  functionality  of  the  groups  that  they 
replace.  We  have  developed  a  model  of  DCPD  molecules  following  this  scheme  and 
ascertained  that  the  polymeric  networks  that  form  have  the  same  geometry  as  those 


(a) 

Fig.  1  Simulated  DCPD  molecule,  (a)  with  all  atoms 
represented  explicitly,  and  (b)  with  hydrogen  atoms  treated 
implicitly.  This  reduced  form  of  the  molecule  contains 
three  types  of  hydrocarbon  croups  sp3  CH2  (red),  sp3  CH 
(black),  and  sp2  CH  (blue).  Reactive  sites  are  the  double 
bonds  between  so2  arouos. 
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resulting  from  all-atom  simulations.  Fig.  2  shows  an  example  of  a  simulation  result  using 
this  simplified  model  of  DCPD  monomers. 


This  simplified  representation  of  the 
DCPD  molecule  results  in  an  in¬ 
crease  of  computational  speed  by 
nearly  an  order  of  magnitude.  Ac¬ 
cordingly,  it  allows  for  an  increase 
in  the  size  of  simulated  structures, 
yielding  improved  statistics  in  the 
description  of  the  resulting  network 
topologies.  Based  on  this  first 
coarse-graining  step,  we  can  gener¬ 
ate  polymeric  networks  ultimately 
containing  several  ten  thousand 
monomers  and  study  their  structure 
and  properties  under  conditions  that 
represent  all  possible  states  the 
healing  agent  can  achieve  between 
initiation  and  completion  of  the  po¬ 
lymerization,  assuming  various  de¬ 
grees  of  mechanical  deformation. 
The  important  measure  that  results 
from  this  level  of  simulation  is  the 
relationship  between  structure  and  mechanical  properties  of  the  polymer  network.  In 
particular,  we  can  determine  elastic  moduli  and  tensile  strengths  of  networked  structures 
at  various  densities.  Variable  density,  in  this  case,  corresponds  to  situations  where 
completion  of  the  polymerization  reaction  was  achieved  at  different  degrees  of  strain  in 
the  system. 


Fig.  2  Fragment  of  a  polymerized  DCPD  network 
as  generated  by  MD  simulations  based  on  the 
unified  building  blocj  approach.  Unreacted 
monomers  are  not  shown  for  claritv. 


3.2  Progress  in  Coarse-Grained  Simulations 

The  simulation  times  that  can  be  achieved  at  the  aforementioned  level,  however,  are  far 
shorter  than  would  be  required  to  cover  the  experimental  deformation  rates.  To 
realistically  accommodate  simultaneous  deformation  and  polymerization  as  the  factors 
contributing  to  the  structural  evolution  of  DCPD  networks,  a  second  level  of  coarse 
graining  has  been  introduced,  i.e.,  the  CGP  procedure.  At  this  level,  each  DCPD 
monomer  is  represented  by  a  single  particle  that  can  move  about,  and  upon  collision  with 
other  units  it  can  undergo  a  polymerization  reaction  depending  on  whether  the  conditions 
for  this  process  to  occur  are  fulfilled.  These  conditions  are  known  from  atomistic 
simulations,  and  accordingly  we  developed  statistical  structural  evolution  rules  that  are 
implemented  via  a  Monte  Carlo  numerical  procedure.  These  rules,  obviously,  are 
designed  to  maintain  realism  in  the  simulated  polymerization  process.  E.g.,  at  least  one 
of  the  particles  must  carry  the  Ru  catalytic  complex.  Upon  reaction,  the  Ru  complex 
residing  with  the  most  recently  polymerized  unit  is  now  transferred  to  the  next  monomer 
that  attaches  to  the  chain.  This  provides  for  an  interesting  aggregation  scheme,  described 
by  the  propagation  of  the  catalyst  along  the  growth  path  of  the  polymer.  Furthermore,  at 


most  four  other  units  can  attach  to  any  given  monomer,  and  we  distinguish  between 
polymerization  and  cross-linking  reactions  based  on  their  respective  rates. 

Initially,  we  tested  this  concept  using  2-D  lattice  Monte  Carlo  model.  This  proof  of 
concept  study  already  yielded  important  design  information.  Specifically,  we  found  that 
for  optimal  utilization  of  the  healing  agent,  it  is  imperative  to  embed  the  catalyst  in  such  a 
way  that  it  dissolves  rapidly  and  be  evenly  distributed  throughout  the  DCPD  fluid  before 
polymerization  reaches  significant  degrees.  Otherwise,  the  catalyst  may  become  trapped 
and  rendered  ineffective.  Alternatively,  it  may  be  necessary  to  slow  the  polymerization 
rate,  e.g.,  by  increasing  the  fraction  of  endo-modification  of  DCPD  in  the  healing  agent. 

We  have  now  completed  the  CGP  method  for  simulation  in  3-D  and  for  particles 
positioned  off  lattice.  Unreacted  monomers  move  according  to  dissipative  particle 
dynamics  algorithm,  which  is  similar  to  an  MD  algorithm,  except  that  interactions 
between  particles  are  designed  to  allow  for  larger  time  steps  at  the  cost  of  reproducing 
only  approximate  particle  trajectories.  This  procedure  is  justified  when  particles 
themselves  are  only  approximations  of  molecules  they  represent.  The  chief  advantage  of 
this  method  is  that  due  to  the  lack  of  conservation  of  dynamics,  the  time  scale  is 
undefined  and  can  be  synchronized  with  those  of  the  reaction  and  deformation  rates  of 
the  experiments.  Synchronization  is  achieved  by  including  a  Monte  Carlo  procedure  that 
accounts  for  polymerization  events.  I.e.,  at  each  particle  collision  a  random  selection 
process  is  invoked  that  results  in  bonding  and  passing  of  the  catalyst  with  a  probability 
that  matches  the  reaction  rate,  as  derived  from  all-atoms  simulations  (or  experiments,  if 
available).  This  defines  the  time  scale,  and  deformation  rates  are  adjusted  accordingly. 


(a)  (b)  (c) 


Fig.  3  (a)  and  (b):  DCPD  structures  at  different  stages  of  the  simulated  polymerization  reaction 
as  evolved  using  our  CGP  method  simulation  of  the  DCPD  polymerization  process.  Grey 
particles  are  unreacted  monomers,  blue  particles  are  polymerized  monomers,  and  green 
particles  are  monomers  carrying  the  catalyst,  (c)  Degree  of  cure  as  a  function  of  time. 

Fig.  3  shows  CGP  simulation  results  for  polymerization  without  deformation.  The  degree 
of  polymerization  changes  with  reaction  time  approximately  like  a  first  order  process. 
The  saturation  level  is  reached  between  85%  and  90%,  depending  on  the  concentration  of 
catalyst.  We  have  furthermore  simulated  this  process  when  subjecting  the  system  to 
various  deformation  rates,  which  directly  yields  the  cohesive  behavior  of  the  healing 


agent,  i.e.,  the  load-deformation  relationship  for  various  cure  rates,  as  needed  as  input  to 
the  CVFE  calculations.  Fig.  4  shows  a  typical  result  for  a  simulation  in  which  the  system 
that  was  cured  to  about  35%  (i.e.,  just  above  the  percolation  threshold)  prior  to  applying 
strain,  and  which  is  allowed  to  contract  laterally  when  undergoing  expansion  in  one 
direction.  Note  that  the  straining  process  begins  at  3.1-1 0-6  s.  The  system  is  already 
under  tensile  stress  at  this  point,  because  it  was  maintained  at  constant  volume  during  the 
initial  curing  phase.  This  initial  stress  can  be  considered  as  a  residual  stress  resulting 
from  the  curing  process.  Fracture  occurs  at  about  230%  strain.  Note  the  necking  near  the 
bottom  of  the  inset  showing  the  structures.  Another  interesting  observation  is  that  the 
reaction  rate  increases  as  soon  as  the  strain  is  applied.  This  is  because  of  the  structural 
rearrangements  during  the  deformation  of  the  polymer  network  bring  unreacted 
monomers  closer  to  one  another. 


t(s) 


Fig.  4  Cohesive  law  of  a  DCPD  system,  subject  to  a  constant  strain  rate  of  10®  s'1.  Insets 
show  the  structures  at  particular  strain  points. 


4.  Direction  of  Future  Work 

Currently  we  are  systematically  investigating  this  process  to  better  understand  the 
statistical  mechanics  of  the  combined  cure  and  deformation  process  in  self-healing 
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polymer  matrix  composites.  We  are  exploring  the  effect  of  catalyst  concentration  and 
mobility,  the  ratio  between  polymerization  and  deformation  rates,  and  the  degree  of  cure 
prior  to  deformation. 

Based  on  these  results  we  will  establish  an  effective  data  connection  with  the  CVFE 
portion  of  this  numerical  framework.  We  will  assess  the  frequency  and  resolution  in 
parameter  space  with  which  the  cohesive  law  for  the  self-healing  polymer  needs  to  be 
communicated  to  the  CVFE  calculations.  We  will  furthermore  develop  an  efficient 
procedure  for  integrating  the  calculations  at  different  scales  into  a  seamless  suite  of 
numerical  techniques. 

Finally,  we  plan  to  develop  a  reverse  procedure  to  the  coarse-graining  at  the  molecular 
scale.  As  outlined  above,  in  the  CGP  simulations  the  description  of  structural  units  and 
the  interactions  between  them  is  simplified.  To  ascertain  that  the  structural  developments 
simulated  by  the  CGP  method  are  sufficiently  accurate,  we  will  develop  a  scheme  for 
mapping  atomically  resolved  structures  onto  the  coarse-grained  ones  in  strategic 
locations,  and  verify  whether  the  coarse-grained  structure  is  compatible  with  atomic 
interactions.  In  case  discrepancies  are  detected,  we  use  this  information  to  iteratively 
correct  the  interaction  and  reaction  rules  at  the  coarse-grained  level. 
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5.  Interactions/Transitions: 


a.  Participation 

AFOSR-MEANS  Bidder’s  Meeting,  May  3, 2004,  Arlington,  VA 

2004  AFOSR  Program  Review  on  “Mechanics  of  Materials  and  Devices.”  August  18-20 
2004,  Wintergreen  Resort,  VA.  Presentation  title:  “Multiscale  Modeling  for  the  Design 
of  Autonomic  Healing  Structural  Composite  Materials.” 


b.  Transitions 

This  work  constitutes  a  part  of  a  multi-investigator  effort  to  bridge  length  and  time  scales 
in  materials  simulation.  Our  group  works  on  the  atomic-scale  and  mesoscale  simulation 
of  the  curing  process  in  the  healing  agent  of  a  self-healing  polymer  composite  material 
under  various  loading  conditions.  The  outcomes  of  our  simulations  are  the  constitutive 
law  for  the  polymer  and  healing  rate  coefficients,  which  serve  as  input  to  the  continuum 
mechanical  calculations  carried  out  by  Prof.  P.  Geubelle  and  his  group  at  the  University 
of  Illinois. 

New  collaborations  with  Dr.  Ajit  Roy  from  the  Air  Force  Research  Laboratory  at  Wright 
Patterson  Air  Force  Base  and  Prof.  Somnath  Gosh  from  Ohio  State  University  have  been 
developed. 
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